package Ascore;

require Exporter;

our $VERSION = '1.00';
our(@ISA, @EXPORT);
@ISA = qw(Exporter);
@EXPORT = qw(ascore conversion_styuv);

use strict;
use warnings;
use Math::Combinatorics;
use Math::NumberCruncher;
use Data::Dumper;


sub ascore{
	my ($peptide,$parent_mass,$z,$masses,$intensities,$precision,$limit_mz_sup,$limit_mz_inf)=@_;
	#print "peptide standard notation: $peptide\n parent: $parent_mass \n charge: $z \n $masses \n $intensities \n";
	#print $peptide,"\n";
	my %result;
	#El numero de sitios de fosforilacion ha de ser mayor que el de 
	#fosforilaciones para que el Ascore tenga sentido. Si el Ascore no tiene sentido, lo devuelvo como 1000
	my $num_sitios_fosforilacion=0;
	my $num_sitios_deshidratacion=0;
	my $num_sitios_modificacion=0;
	my $num_fosforilaciones=0;
	my $num_deshidrataciones=0;
	my $num_modificaciones=0;
	
	my $peptido_conteo_modificaciones=$peptide;
	while($peptido_conteo_modificaciones=~s/(\(21\)|\(21\,|21\))/x/){
		$num_fosforilaciones++;
	}
	while($peptido_conteo_modificaciones=~s/(\(23\)|\(23\,|23\))/x/){
		$num_deshidrataciones++;
	}
	while($peptido_conteo_modificaciones=~s/^Y/x/){
		$num_sitios_fosforilacion++;
		$num_sitios_modificacion++;
	}
	while($peptido_conteo_modificaciones=~s/^[ST]/x/){
		$num_sitios_fosforilacion++;
		$num_sitios_deshidratacion++;
		$num_sitios_modificacion++;
	}
	while($peptido_conteo_modificaciones=~s/Y(\w)/x$1/){
		$num_sitios_fosforilacion++;
		$num_sitios_modificacion++;
	}
	while($peptido_conteo_modificaciones=~s/[ST](\w)/x$1/){
		$num_sitios_fosforilacion++;
		$num_sitios_deshidratacion++;
		$num_sitios_modificacion++;
	}

	#Si el numero de sitios de modificacion es el mismo que el de modificaciones
	$num_modificaciones = $num_fosforilaciones + $num_deshidrataciones;
	if($num_sitios_modificacion==$num_modificaciones){
		$result{'phosphorylations'}=$num_fosforilaciones;
		$result{'dehydratations'}=$num_deshidrataciones;
		$result{'peptide'}=$peptide;
		if($num_modificaciones==1){
			$result{'ascore'}=[1000];
		}
		elsif($num_modificaciones==2){
			$result{'ascore'}=[1000,999];
		}
		elsif($num_modificaciones==3){
			$result{'ascore'}=[1000,999,998];
		}
		return %result;
	}
	#Caso especial: numero de deshidrataciones es igual al numero de posibles sitios de deshidratacion
		if($num_deshidrataciones>0){
			if($num_sitios_deshidratacion==$num_deshidrataciones){
				if($num_modificaciones==1){
					$result{'phosphorylations'}=$num_fosforilaciones;
					$result{'dehydratations'}=$num_deshidrataciones;
					$result{'peptide'}=$peptide;
					$result{'ascore'}=[1000];
					return %result;
				}
				else{
					#Caso no resuelto: num.DHs == num.sitios.DH Y num_mods>=2
					$result{'phosphorylations'}=$num_fosforilaciones;
					$result{'dehydratations'}=$num_deshidrataciones;
					$result{'peptide'}=$peptide;
					if($num_modificaciones==2){
						$result{'ascore'}=[0,0];
					}
					elsif($num_modificaciones==3){
						$result{'ascore'}=[0,0,0];
					}
					print "CASO NO RESUELTO \n";
				}
				
			}
		}
	if($num_modificaciones==0){
		$result{'phosphorylations'} = $num_fosforilaciones;
		$result{'dehydratations'} = $num_deshidrataciones;
		$result{'peptide'} = $peptide;
	}
	elsif($num_modificaciones==1){
		%result = ascore_1_fosforilacion($peptide,$precision,$masses,$intensities,$z,$limit_mz_sup,$limit_mz_inf);
	}
	elsif($num_modificaciones==2){
		%result = ascore_2_fosforilaciones($peptide,$precision,$masses,$intensities,$z,$limit_mz_sup,$limit_mz_inf);
	}
	elsif($num_modificaciones==3){
		%result = ascore_3_fosforilaciones($peptide,$precision,$masses,$intensities,$z,$limit_mz_sup,$limit_mz_inf);
	}
	else{
		$result{'phosphorylations'}=$num_fosforilaciones;
		$result{'dehydratations'}=$num_deshidrataciones;
		$result{'peptide'}=$peptide;
		
	}
	return %result;
}

sub ascore_1_fosforilacion{
	my $peptido = $_[0];			#en formato [styuv]
	my $precision = $_[1];			#en unidades m/z
	my $masses = $_[2];				#referencia a array de masas
	my $intensities = $_[3];		#referencia a array intensidades
	my $z = $_[4];					#carga
	my $limit_mz_sup = $_[5];			
	my $limit_mz_inf = $_[6];

	my @masas = split /\|/,$masses;
	my @intensidades = split /\|/,$intensities;
	
	#print $peptido,"\n";#exit;
	my @peptidos_posibles = peptide_collection($peptido);
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z,$limit_mz_sup,$limit_mz_inf)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score = $peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{'score'}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{'score_ponderado'} = ($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{'score'}/$suma_niveles);
		}
		$peptide_scores{$pep} = $peptide_score;
	}
	
	#Ordeno peptidos segun Peptide score
	my @peptidos_con_sumas_scores_ordenadas=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$peptidos_con_sumas_scores_ordenadas[0];
	my $p2=$peptidos_con_sumas_scores_ordenadas[1];
	my $p3=$peptidos_con_sumas_scores_ordenadas[2];
	
	#Convierto p1 y p2 a peptidos en formato [styuv]
	my $p1_s=conversion_styuv($p1);
	my $p2_s=conversion_styuv($p2);

	#print $p1_s,"\n",$p2_s,"\n";
	#if(($peptide_scores{$p2})and($peptide_scores{$p1}==$peptide_scores{$p2})){
		#print "Peptide scores iguales para peptidos 1 y 2\n";
		#return ($peptido,$peptide_scores{$p1},0,0);
	#}
	#testeo si, en el caso de haber un tercer candidato, su peptide score es tan bueno como el segundo;
	if (defined $p3){ 
		my $p3_s=conversion_styuv($p3);			
		if((defined $peptide_scores{$p3})and($peptide_scores{$p2}==$peptide_scores{$p3})){
			#Si p2 tiene el mismo Peptide Score que p3, debo elegor uno para usar en el ascore 
			#print "pep2 ($p2) y pep 3 ($p3) son iguales ($peptide_scores{$p2}), respecto a $p1\n";
			###Determino cual de los peptidos p2 y p3 tiene la fosforilacion mas proxima a p1
			my @p_trozos_peptido_1= split //,$p1_s;
			my @p_trozos_peptido_2= split //,$p2_s;
			my @p_trozos_peptido_3= split //,$p3_s;
			my ($punto1,$punto2,$punto3);
			for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){if($p_trozos_peptido_1[$i]=~/[styuv]/){$punto1=1+$i;}}
			for(my $i=0;$i<$#p_trozos_peptido_2+1;$i++){if($p_trozos_peptido_2[$i]=~/[styuv]/){$punto2=1+$i;}}
			for(my $i=0;$i<$#p_trozos_peptido_3+1;$i++){if($p_trozos_peptido_3[$i]=~/[styuv]/){$punto3=1+$i;}}
			my $p2_usar; #uso el peptido con menor diferencia a punto1
			$punto2=abs($punto1-$punto2);
			$punto3=abs($punto1-$punto3);
			if($punto2<$punto3){
				$p2_usar=$p2;
			}
			else{
				$p2_usar=$p3;
			}
			$p2=$p2_usar;
			#print "Peptido usado: ".$p2_usar."\n";
			
		}
	}
	
	my @peptidos_ascore;
	push @peptidos_ascore,$p1;
	push @peptidos_ascore,$p2;
	
	#Obtengo nivel maxima diferencia entre Peptide Scores
	my $level_max_diff=0;
	my $difss=0;
	for(my $i=1;$i<11;$i++){
		my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$p2}{$i}{score};
		if($dif>$difss){
			$difss=$dif;
			$level_max_diff=$i;
		}
	}
	if($level_max_diff==0){
		$level_max_diff=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
	}
	
	my @p_trozos_peptido_1= split //,$p1_s;
	my @p_trozos_peptido_2= split //,$p2_s;
	my ($punto1,$punto2);
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			$punto1=1+$i;
		}
	}
	for(my $i=0;$i<$#p_trozos_peptido_2+1;$i++){
		if($p_trozos_peptido_2[$i]=~/[styuv]/){
			$punto2=1+$i;
		}
	}
	if($punto1>$punto2){
		($punto1,$punto2)=($punto2,$punto1);
	}
	
	my %fragmentos_p1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p1=%{tratamiento_fragmentos($p1_s,\%fragmentos_p1,$punto1,$punto2)};
	my %fragmentos_p2=%{prever_fragmentos($p2,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p2=%{tratamiento_fragmentos($p2_s,\%fragmentos_p2,$punto1,$punto2)};

	
	my %datos_ascore;
	$datos_ascore{$p1}=matches_ascore(\%fragmentos_p1,\@masas,\@intensidades,$precision,$p1);
	$datos_ascore{$p2}=matches_ascore(\%fragmentos_p2,\@masas,\@intensidades,$precision,$p2);
			
	my $ascore;
	#print "p1-> $datos_ascore{$p1}{$level_max_diff}{score}\n";
	#print "p2-> $datos_ascore{$p2}{$level_max_diff}{score}\n";
	
	my $matches1=1+$#{$datos_ascore{$p1}{$level_max_diff}{tags_match_intervalo}};
	my $matches2=1+$#{$datos_ascore{$p2}{$level_max_diff}{tags_match_intervalo}};
	my @array1=keys(%fragmentos_p1);
	my $size=1+$#array1;
	my $proporcion=$matches1."//".$size." (".$matches2.")";
	
	#if((not defined $datos_ascore{$p1}{$level_max_diff}{score})and(not defined $datos_ascore{$p2}{$level_max_diff}{score})){
		#print "No definido score p1 ni score p2\n";
		#return ($p1,$level_max_diff,$peptide_scores{$p1},0,$proporcion);
	#}
	#if((not defined $datos_ascore{$p1}{$level_max_diff}{score})){
		#print "No definido score p1\n";
		#return ($p1,$level_max_diff,$peptide_scores{$p1},0,$proporcion);
	#}

	$ascore=($datos_ascore{$p1}{$level_max_diff}{score})-($datos_ascore{$p2}{$level_max_diff}{score});	
	#return ($p1,$level_max_diff,$peptide_scores{$p1},$p2."(".$punto1."-".$punto2.")"," ",$ascore,$proporcion);
	my %result;
	my $num_deshidrataciones=0;
	my $pdh=$p1;
		while($pdh=~s/23//){
			$num_deshidrataciones++;
		}
	my $num_fosforilaciones=0;
	my $pph=$p1;
		while($pph=~s/(\(21\)|\(21\,|21\))//){
			$num_fosforilaciones++;
		}
	$result{'phosphorylations'}=$num_fosforilaciones;
	$result{'dehydratations'}=$num_deshidrataciones;
	$result{'peptide'}=$p1;
	$result{'peptide_score'}=$peptide_scores{$p1};
	if($ascore<0.1){
		$ascore=0;
	}
	$result{'ascore'}=[$ascore];
	
	return %result;

}


sub ascore_2_fosforilaciones{
	#print "2 fosfos\n";
	my $peptido = $_[0];			#en formato [styuv]
	my $precision = $_[1];		#en unidades m/z
	my $masses = $_[2];			#referencia a array de masas
	my $intensities = $_[3];	#referencia a array intensidades
	my $z = $_[4];				#carga
	my $limit_mz_sup = $_[5];			#Limite del maximo valor mz usado (2000?)
	my $limit_mz_inf = $_[6];
	#print "Ascore 2 fosfos\n";
	my @masas = split /\|/,$masses;
	my @intensidades = split /\|/,$intensities;

	
	my @peptidos_posibles = peptide_collection($peptido); 
	
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos= %{prever_fragmentos($pep,$z,$limit_mz_sup,$limit_mz_inf)};
		$datos_matches_10_niveles{$pep}=matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	

	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score=$peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{score_ponderado}=($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}

	#Ordeno peptidos segun Peptide score
	my @pepts_sorted_by_score=sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1=$pepts_sorted_by_score[0];
	
	#Convierto peptidos en formato [styuv] segun se vaya necesitando
	
	#Detecto puntos de modificacion
	my @p_trozos_peptido_1= split //,conversion_styuv($p1);
	my $f1=0;
	my $f2=0;
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			my $punto=1+$i;
			if($f1){
				$f2=$punto;
			}
			else{
				$f1=$punto;
			}
		}
	}
	
	#Encuentro Peptidos-2 para cada modificacion
	my $peptido2_f1;
	my $peptido2_f2;
	for (my $i=1;$i<$#pepts_sorted_by_score+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($pepts_sorted_by_score[$i]);
		if(($p_trozos_peptido_x[(-1+$f1)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)){
			$peptido2_f1=$pepts_sorted_by_score[$i];
			last;
		}
	}
	for (my $i=1;$i<$#pepts_sorted_by_score+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($pepts_sorted_by_score[$i]);
		if(($p_trozos_peptido_x[(-1+$f2)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)){
			$peptido2_f2=$pepts_sorted_by_score[$i];
			last;
		}
	}

	#Obtengo niveles maxima diferencia entre Peptide Scores (P1 respecto a peptido2_f1 y peptido2_f2)
	my $level_max_diff1=0;
	my $difss1=0;
	my $level_max_diff2=0;
	my $difss2=0;
	{
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f1}{$i}{score};
			if($dif>$difss1){
				$difss1=$dif;
				$level_max_diff1=$i;
			}
		}
		if($level_max_diff1==0){
			$level_max_diff1=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f2}{$i}{score};
			if($dif>$difss2){
				$difss2=$dif;
				$level_max_diff2=$i;
			}
		}
		if($level_max_diff2==0){
			$level_max_diff2=5;#Esto ocurre cuando los scores del pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
	}
	
	#Obtengo fragmentos dentro de los intervalos adecuados
		#1#Encuentro puntos complementarios a las f1 y f2 en cada peptido
		my @p_trozos_peptido_2f1= split //,conversion_styuv($peptido2_f1);
		my @p_trozos_peptido_2f2= split //,conversion_styuv($peptido2_f2);
		my $aa_complement_1;
		for(my $i=0;$i<$#p_trozos_peptido_2f1+1;$i++){
			if($p_trozos_peptido_2f1[$i]=~/[a-z]/){
				unless((1+$i)==$f2){
					$aa_complement_1=1+$i;
				}
			}
		}
		my $aa_complement_2;
		for(my $i=0;$i<$#p_trozos_peptido_2f2+1;$i++){
			if($p_trozos_peptido_2f2[$i]=~/[a-z]/){
				unless((1+$i)==$f1){
					$aa_complement_2=1+$i;
				}
			}
		}
		
		#2#Ordeno los puntos que marcan los intervalos adecuados
		my ($punto11,$punto12,$punto21,$punto22);
		if($f1>$aa_complement_1){
			$punto11=$aa_complement_1;
			$punto12=$f1;
		}
		elsif($f1<$aa_complement_1){
			$punto12=$aa_complement_1;
			$punto11=$f1;
		}
		else{
			print "errorrrrrrrrrrrrrrr!!!!\n";
		}
		
		if($f2>$aa_complement_2){
			$punto21=$aa_complement_2;
			$punto22=$f2;
		}
		elsif($f2<$aa_complement_2){
			$punto22=$aa_complement_2;
			$punto21=$f2;
		}
		else{
			print "errorrrrrrrrrrrrrrr!!!!\n";
		}
		
		
		#3#Calculo fragmentos
		#Fragmentos para fosforilacion 1
		my %fragmentos_p1_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
		%fragmentos_p1_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p1_1,$punto11,$punto12)};
		my %fragmentos_p1_2=%{prever_fragmentos($peptido2_f1,$z,$limit_mz_sup,$limit_mz_inf)};
		%fragmentos_p1_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f1),\%fragmentos_p1_2,$punto11,$punto12)};
		#Fragmentos para fosforilacion 2
		my %fragmentos_p2_1=%{prever_fragmentos($p1,$z,$limit_mz_sup,$limit_mz_inf)};
		%fragmentos_p2_1=%{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p2_1,$punto21,$punto22)};
		my %fragmentos_p2_2=%{prever_fragmentos($peptido2_f2,$z,$limit_mz_sup,$limit_mz_inf)};
		%fragmentos_p2_2=%{tratamiento_fragmentos(conversion_styuv($peptido2_f2),\%fragmentos_p2_2,$punto21,$punto22)};
		
		##Obtengo datos-ascore: con matches de los fragmentos site-determining
		my $datos_ascore_p1_1=matches_ascore(\%fragmentos_p1_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p1_2=matches_ascore(\%fragmentos_p1_2,\@masas,\@intensidades,$precision,$peptido2_f1);
				
		my $datos_ascore_p2_1=matches_ascore(\%fragmentos_p2_1,\@masas,\@intensidades,$precision,$p1);
		my $datos_ascore_p2_2=matches_ascore(\%fragmentos_p2_2,\@masas,\@intensidades,$precision,$peptido2_f2);
		
		my $matches1_1=1+$#{${$datos_ascore_p1_1}{$level_max_diff1}{tags_match_intervalo}};
		my $matches1_2=1+$#{${$datos_ascore_p1_2}{$level_max_diff1}{tags_match_intervalo}};
		my @array1=keys(%fragmentos_p1_1);
		my $size1=1+$#array1;
		my $proporcion1=$matches1_1."//".$size1." (".$matches1_2.")";
		
		my $matches2_1=1+$#{${$datos_ascore_p2_1}{$level_max_diff2}{tags_match_intervalo}};
		my $matches2_2=1+$#{${$datos_ascore_p2_2}{$level_max_diff2}{tags_match_intervalo}};
		my @array2=keys(%fragmentos_p2_1);
		my $size2=1+$#array2;
		my $proporcion2=$matches2_1."//".$size2." (".$matches2_2.")";
	
		my $ascore1 = ${$datos_ascore_p1_1}{$level_max_diff1}{score}-${$datos_ascore_p1_2}{$level_max_diff1}{score};	
		my $ascore2 = ${$datos_ascore_p2_1}{$level_max_diff2}{score}-${$datos_ascore_p2_2}{$level_max_diff2}{score};
		
		
		#print $p1," -> ";
		#print $ascore1,"\t|\t",$ascore2,"\n";
		#print "Level_max_1  $level_max_diff1\n";
		#print "Level_max_2  $level_max_diff2\n";
		#print "peptido de referencia  $p1\n";
		#print "pept1 = $peptido2_f1\n";
		#print "pept2 = $peptido2_f2\n";
		
		
	my %result;
	my $num_deshidrataciones=0;
	my $pdh=$p1;
		while($pdh=~s/23//){
			$num_deshidrataciones++;
		}
	my $num_fosforilaciones=0;
	my $pph=$p1;
		while($pph=~s/(\(21\)|\(21\,|21\))//){
			$num_fosforilaciones++;
		}
		
	$result{'phosphorylations'}=$num_fosforilaciones;
	$result{'dehydratations'}=$num_deshidrataciones;
	$result{'peptide'}=$p1;
	$result{'peptide_score'}=$peptide_scores{$p1};
	if($ascore1<0.1){
		$ascore1 = 0;
	}
	if($ascore2<0.1){
		$ascore2 = 0;
	}
	$result{'ascore'} = [$ascore1, $ascore2];
	
	return %result;	
}


sub ascore_3_fosforilaciones{
	my $peptido = $_[0];		#en formato [styuv]
	my $precision = $_[1];		#en unidades m/z
	my $masses = $_[2];			#referencia a array de masas
	my $intensities = $_[3];	#referencia a array intensidades
	my $z = $_[4];				#carga
	my $limit_mz_sup = $_[5];			
	my $limit_mz_inf = $_[6];
	#print "Ascore 2 fosfos\n";
	my @masas = split /\|/,$masses;
	my @intensidades = split /\|/,$intensities;

	my @peptidos_posibles = peptide_collection($peptido); 
	
	#Obtencion de los Peptide Scores
	my %peptide_scores;
	my %datos_matches_10_niveles;
	foreach my $pep(@peptidos_posibles){
		my %fragmentos_teoricos = %{prever_fragmentos($pep, $z, $limit_mz_sup, $limit_mz_inf)};
		$datos_matches_10_niveles{$pep} = matches_10niveles_series_yb_completas(\%fragmentos_teoricos,\@masas,\@intensidades,$precision,$pep);
	}
	
	my %weigths=(1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25,	10=>0.25);
	my $suma_niveles=7;	#suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25
	

	foreach my $pep(keys %datos_matches_10_niveles){
		my $suma_datos_matches=0;
		my $peptide_score=0;
		for (my $i=1;$i<11;$i++){
			$peptide_score = $peptide_score+($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
			$datos_matches_10_niveles{$pep}{$i}{score_ponderado} = ($weigths{$i}*$datos_matches_10_niveles{$pep}{$i}{score}/$suma_niveles);
		}
		$peptide_scores{$pep}=$peptide_score;
	}

	#Ordeno peptidos segun Peptide score
	my @pepts_sorted_by_score = sort {$peptide_scores{$b} <=> $peptide_scores{$a}} keys %peptide_scores;
	my $p1 = $pepts_sorted_by_score[0];
	
	#Convierto peptidos en formato [styuv] segun se vaya necesitando
	
	#Detecto puntos de modificacion
	my @p_trozos_peptido_1= split //,conversion_styuv($p1);
	my $f1=0;
	my $f2=0;
	my $f3=0;
	for(my $i=0;$i<$#p_trozos_peptido_1+1;$i++){
		if($p_trozos_peptido_1[$i]=~/[styuv]/){
			my $punto=1+$i;
			if($f2){
				$f3=$punto;
			}
			elsif($f1){
				$f2=$punto;
			}
			else{
				$f1=$punto;
			}
		}
	}

	#Encuentro Peptidos-2 para cada modificacion
	#Los peptidos-2 se escogen con este criterio: se busca un peptido con el mayor peptide-score posible y que
	#	cumpla la condicion que el p1 NO este modificado y que los otros p1(con 3 modificaciones, para el caso de
	# 	buscar para f1, el p1 de f2 y el p1 de f3) SI esten modificados.	
	my $peptido2_f1;
	my $peptido2_f2;
	my $peptido2_f3;
	for (my $i=1;$i<$#pepts_sorted_by_score+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($pepts_sorted_by_score[$i]);
		if(($p_trozos_peptido_x[(-1+$f1)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)and($p_trozos_peptido_x[(-1+$f3)]=~/[a-z]/)){
			$peptido2_f1=$pepts_sorted_by_score[$i];
			last;
		}
	}
	for (my $i=1;$i<$#pepts_sorted_by_score+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($pepts_sorted_by_score[$i]);
		if(($p_trozos_peptido_x[(-1+$f2)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)and ($p_trozos_peptido_x[(-1+$f3)]=~/[a-z]/)){
			$peptido2_f2=$pepts_sorted_by_score[$i];
			last;
		}
	}
	for (my $i=1;$i<$#pepts_sorted_by_score+1;$i++){
		my @p_trozos_peptido_x=split //,conversion_styuv($pepts_sorted_by_score[$i]);
		if(($p_trozos_peptido_x[(-1+$f3)]=~/[A-Z]/)and($p_trozos_peptido_x[(-1+$f1)]=~/[a-z]/)and ($p_trozos_peptido_x[(-1+$f2)]=~/[a-z]/)){
			$peptido2_f3=$pepts_sorted_by_score[$i];
			last;
		}
	}

	#Obtengo niveles maxima diferencia entre Peptide Scores (P1 respecto a peptido2_f1, peptido2_f2 y peptido2_f3)
	my $level_max_diff1 = 0;
	my $difss1 = 0;
	my $level_max_diff2=0;
	my $difss2 = 0;
	my $level_max_diff3=0;
	my $difss3 = 0;
	{
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f1}{$i}{score};
			if($dif>$difss1){
				$difss1=$dif;
				$level_max_diff1=$i;
			}
		}
		if($level_max_diff1==0){
			$level_max_diff1=5;#Esto ocurre cuando el pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f2}{$i}{score};
			if($dif>$difss2){
				$difss2=$dif;
				$level_max_diff2=$i;
			}
		}
		if($level_max_diff2==0){
			$level_max_diff2=5;#Esto ocurre cuando los scores del pep1 y el pep2 son identicos; uso el 5 de forma arbitraria 
		}
		for(my $i=1;$i<11;$i++){
			my $dif= $datos_matches_10_niveles{$p1}{$i}{score}-$datos_matches_10_niveles{$peptido2_f3}{$i}{score};
			if($dif>$difss3){
				$difss3=$dif;
				$level_max_diff1=$i;
			}
		}
		if($level_max_diff3==0){
			$level_max_diff3=5;#Esto ocurre cuando el pep1 y el pep3 son identicos; uso el 5 de forma arbitraria 
		}
	}
	
	#Obtengo fragmentos dentro de los intervalos adecuados
	#   1#Encuentro puntos complementarios a las f1,f2 y f3 en cada peptido
	#	Por punto complementario se entiende al aminoacido que este modificado en p2 pero NO en p1
	my @p_trozos_peptido_2f1 = split //,conversion_styuv($peptido2_f1);
	my @p_trozos_peptido_2f2 = split //,conversion_styuv($peptido2_f2);
	my @p_trozos_peptido_2f3 = split //,conversion_styuv($peptido2_f3);
	
	my $aa_complement_1;
	for(my $i=0; $i<$#p_trozos_peptido_2f1+1; $i++){
		if($p_trozos_peptido_2f1[$i]=~/[a-z]/){
			unless(((1+$i)==$f2)or((1+$i)==$f3)){
				$aa_complement_1 = 1+$i;
			}
		}
	}
	my $aa_complement_2;
	for(my $i=0;$i<$#p_trozos_peptido_2f2+1;$i++){
		if($p_trozos_peptido_2f2[$i]=~/[a-z]/){
			unless(((1+$i)==$f1)or((1+$i)==$f3)){
				$aa_complement_2=1+$i;
			}
		}
	}
	my $aa_complement_3;
	for(my $i=0;$i<$#p_trozos_peptido_2f3+1;$i++){
		if($p_trozos_peptido_2f3[$i]=~/[a-z]/){
			unless(((1+$i)==$f1)or((1+$i)==$f2)){
				$aa_complement_3=1+$i;
			}
		}
	}
	
	#   2#Ordeno los puntos que marcan los intervalos adecuados
	#	Esto significa que ordeno los pares f1-aa_complement1, pares f2-aa_complement2 y 
	#	pares f3-aa_complement3 como pares menor-mayor (sea cual sea el de p1 y p2)
	   
	my ($punto11, $punto12, $punto21, $punto22, $punto31, $punto32);
	if($f1>$aa_complement_1){
		$punto11 = $aa_complement_1;
		$punto12 = $f1;
	}
	elsif($f1<$aa_complement_1){
		$punto12 = $aa_complement_1;
		$punto11 = $f1;
	}
	else{
		print "errorrrrrrrrrrrrrrr!!!!\n";
	}
	
	if($f2>$aa_complement_2){
		$punto21 = $aa_complement_2;
		$punto22 = $f2;
	}
	elsif($f2<$aa_complement_2){
		$punto22 = $aa_complement_2;
		$punto21 = $f2;
	}
	else{
		print "errorrrrrrrrrrrrrrr!!!!\n";
	}
	if($f3>$aa_complement_3){
		$punto31 = $aa_complement_3;
		$punto32 = $f3;
	}
	elsif($f3<$aa_complement_3){
		$punto32 = $aa_complement_3;
		$punto31 = $f3;
	}
	else{
		print "errorrrrrrrrrrrrrrr!!!!\n";
	}
	
	#3#Calculo fragmentos
	#Fragmentos para fosforilacion 1
	my %fragmentos_p1_1 = %{prever_fragmentos($p1, $z, $limit_mz_sup, $limit_mz_inf)};
	%fragmentos_p1_1 = %{tratamiento_fragmentos(conversion_styuv($p1),\%fragmentos_p1_1,$punto11, $punto12)};
	my %fragmentos_p1_2 = %{prever_fragmentos($peptido2_f1,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p1_2 = %{tratamiento_fragmentos(conversion_styuv($peptido2_f1),\%fragmentos_p1_2,$punto11, $punto12)};
	
	#Fragmentos para fosforilacion 2
	my %fragmentos_p2_1 = %{prever_fragmentos($p1, $z, $limit_mz_sup, $limit_mz_inf)};
	%fragmentos_p2_1 = %{tratamiento_fragmentos(conversion_styuv($p1), \%fragmentos_p2_1,$punto21, $punto22)};
	my %fragmentos_p2_2 = %{prever_fragmentos($peptido2_f2,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p2_2 = %{tratamiento_fragmentos(conversion_styuv($peptido2_f2), \%fragmentos_p2_2,$punto21, $punto22)};
	
	#Fragmentos para fosforilacion 3
	my %fragmentos_p3_1 = %{prever_fragmentos($p1, $z, $limit_mz_sup, $limit_mz_inf)};
	%fragmentos_p3_1 = %{tratamiento_fragmentos(conversion_styuv($p1), \%fragmentos_p3_1, $punto31, $punto32)};
	my %fragmentos_p3_2 = %{prever_fragmentos($peptido2_f3,$z,$limit_mz_sup,$limit_mz_inf)};
	%fragmentos_p3_2 = %{tratamiento_fragmentos(conversion_styuv($peptido2_f3),\%fragmentos_p3_2,$punto31,$punto32)};
	
	
	##Obtengo datos-ascore: con matches de los fragmentos site-determining
	my $datos_ascore_p1_1 = matches_ascore(\%fragmentos_p1_1, \@masas, \@intensidades, $precision, $p1);
	my $datos_ascore_p1_2 = matches_ascore(\%fragmentos_p1_2, \@masas, \@intensidades, $precision, $peptido2_f1);
			
	my $datos_ascore_p2_1 = matches_ascore(\%fragmentos_p2_1, \@masas, \@intensidades, $precision, $p1);
	my $datos_ascore_p2_2 = matches_ascore(\%fragmentos_p2_2, \@masas, \@intensidades, $precision, $peptido2_f2);
	
	my $datos_ascore_p3_1 = matches_ascore(\%fragmentos_p3_1, \@masas, \@intensidades, $precision, $p1);
	my $datos_ascore_p3_2 = matches_ascore(\%fragmentos_p3_2, \@masas, \@intensidades, $precision, $peptido2_f3);
	
	#print Dumper $datos_ascore_p1_1;
	
	
	#print $#{${$datos_ascore_p1_1}{$level_max_diff1}{tags_match_intervalo}};
	#exit;
	my $matches1_1=1+$#{${$datos_ascore_p1_1}{$level_max_diff1}{tags_match_intervalo}};
	my $matches1_2=1+$#{${$datos_ascore_p1_2}{$level_max_diff1}{tags_match_intervalo}};
	my @array1=keys(%fragmentos_p1_1);
	my $size1=1+$#array1;
	my $proporcion1=$matches1_1."//".$size1." (".$matches1_2.")";
	
	my $matches2_1 = 1+$#{${$datos_ascore_p2_1}{$level_max_diff2}{tags_match_intervalo}};
	my $matches2_2 = 1+$#{${$datos_ascore_p2_2}{$level_max_diff2}{tags_match_intervalo}};
	my @array2=keys(%fragmentos_p2_1);
	my $size2 = 1+$#array2;
	my $proporcion2=$matches2_1."//".$size2." (".$matches2_2.")";
	
	my $matches3_1 = 1+$#{${$datos_ascore_p3_1}{$level_max_diff3}{tags_match_intervalo}};
	my $matches3_2 = 1+$#{${$datos_ascore_p3_2}{$level_max_diff3}{tags_match_intervalo}};
	my @array3 = keys(%fragmentos_p3_1);
	my $size3 = 1+$#array3;
	my $proporcion3 = $matches3_1."//".$size3." (".$matches3_2.")";
	
	
	my $ascore1 = ${$datos_ascore_p1_1}{$level_max_diff1}{score} -
				  ${$datos_ascore_p1_2}{$level_max_diff1}{score};	
	my $ascore2 = ${$datos_ascore_p2_1}{$level_max_diff2}{score} -
				  ${$datos_ascore_p2_2}{$level_max_diff2}{score};	
	my $ascore3 = ${$datos_ascore_p3_1}{$level_max_diff3}{score} -
				  ${$datos_ascore_p3_2}{$level_max_diff3}{score};
	#print $p1," -> ";
	#print $ascore1,"\t|\t",$ascore2,"\n";
	
	my %result;
	my $num_deshidrataciones=0;
	my $pdh = $p1;
		while($pdh=~s/23//){
			$num_deshidrataciones++;
		}
	my $num_fosforilaciones=0;
	my $pph=$p1;
		while($pph=~s/(\(21\)|\(21\,|21\))//){
			$num_fosforilaciones++;
		}
		
	$result{'phosphorylations'} = $num_fosforilaciones;
	$result{'dehydratations'} = $num_deshidrataciones;
	$result{'peptide'} = $p1;
	$result{'peptide_score'} = $peptide_scores{$p1};
	if($ascore1<0.1){
		$ascore1 = 0;
	}
	if($ascore2<0.1){
		$ascore2 = 0;
	}
	if($ascore3<0.1){
		$ascore3 = 0;
	}
	$result{'ascore'} = [$ascore1, $ascore2, $ascore3];
	
	return %result;	
	
}

sub peptide_collection{
	my $pep = $_[0];
	#Asumo dos cosas:
	#	1/Solo hay un chemtag por cada aminoacido
	#	2/Los chemtag no cambian de sitio
	my @posibles_puntos_fosforilacion;
	my $num_fosforilaciones = 0;
	my @posibles_puntos_DH;
	my $num_DH = 0;
	my $num_chemtag = 0;

	my @aas;
	my @mods;
	my $aa_num = 0;

	while($pep=~/\w/){
		my $aa;
		my $mod;
		if($pep=~s/^(\w)\((\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
			#Calculo el numero de fosfos/chemtags/DH
			if($mod==21){
				$num_fosforilaciones++;
			}
			elsif($mod==23){
				$num_DH++;
			}
			elsif($mod==214 or $mod==737){
				$num_chemtag++;
				if($aa_num>0){#Si hay chemtag en un aa no N-terminal, no se puede modificar con nada mas
					$aa_num++;
					next;
				}
			}
			#Recopilo los posibles puntos de fosfo/DH
			if($aa=~/[S,T,Y]/){
				push @posibles_puntos_fosforilacion,$aa_num;
			}
			if($aa=~/[S,T]/){
				push @posibles_puntos_DH,$aa_num;
			}
		}
		else{
			$pep=~s/^(\w)//;
			$aa=$1;
			$mod='';
			push @aas,$aa;
			push @mods,$mod;
			#Recopilo los posibles puntos de fosfo/chemtag/DH
			if($aa=~/[S,T,Y]/){
				push @posibles_puntos_fosforilacion,$aa_num;
			}
			if($aa=~/[S,T]/){
				push @posibles_puntos_DH,$aa_num;
			}
		}
		$aa_num++;
	}
		
	#Obtengo las diferentes combinaciones de fosforilados, deshidratados y chemtag
	my @combinaciones_fosfo;
	my @combinaciones_DH;
	my @combinaciones_chemtag;
	
	my $combinat1 = Math::Combinatorics->new(count => $num_fosforilaciones,data => [@posibles_puntos_fosforilacion],);
	while(my @combo1 = $combinat1->next_combination){
		push @combinaciones_fosfo,\@combo1;
	}
	my $combinat2 = Math::Combinatorics->new(count => $num_DH,data => [@posibles_puntos_DH],);
		while(my @combo2 = $combinat2->next_combination){
			push @combinaciones_DH,\@combo2;
	}

	my @combinaciones_peptidos;
	if(@combinaciones_fosfo and @combinaciones_DH){
		#print "Hay phospho, DH y chemtag\t";
		foreach my $fos(@combinaciones_fosfo){
			my %usados;
			my @fos=@$fos;
			foreach my $f(@fos){
				$usados{$f}=1;
			}
			LABEL1:
			foreach my $dh(@combinaciones_DH){
				my @dhs=@$dh;
				foreach my $d(@dhs){
					if(defined $usados{$d}){
						next LABEL1;
					}
				}
				#Genero un nuevo array de modificaciones para construir el peptido modificado
				my @mods_new;
				for(my $i=0;$i<$#mods+1;$i++){#Recupero las modificaciones no phospho o DH (chemtag,MetOx,...)
					if($mods[$i]){
						unless(($mods[$i]==21)or($mods[$i]==23)){
							$mods_new[$i]='('.$mods[$i].')';
						}
					}
				}
				foreach(@fos){
					if($mods_new[$_]){
						my $t=$mods_new[$_];
						$t=~s/[\(,\)]//g;
						$mods_new[$_]='('.$t.',21)';	
					}
					else{
						$mods_new[$_]='(21)';
					}
				}
				foreach(@dhs){
					if($mods_new[$_]){
						my $t=$mods_new[$_];
						$t=~s/[\(,\)]//g;
						$mods_new[$_]='('.$t.',23)';	
					}
					else{
						$mods_new[$_]='(23)';
					}
				}
				#Genero el peptido
				my $peptide_string = '';
				for(my $i=0; $i<$#aas+1; $i++){
					if($mods_new[$i]){
						$peptide_string=$peptide_string.$aas[$i].$mods_new[$i];
					}
					else{
						$peptide_string=$peptide_string.$aas[$i];
					}
				}
				push @combinaciones_peptidos,$peptide_string;
			}
		}
	}
	elsif(@combinaciones_fosfo){
		#print "solo fosfo\t";
		foreach my $fos(@combinaciones_fosfo){
			my @fos=@$fos;
			# Nuevo array de modificaciones para construir el peptido modificado
			my @mods_new;
			# Recupero las modificaciones no phospho o DH (chemtag,MetOx,...)
			for(my $i=0;$i<$#mods+1;$i++){
				if($mods[$i]){
					unless(($mods[$i]==21)or($mods[$i]==23)){
						$mods_new[$i] = '('.$mods[$i].')';
					}
				}
			}
			foreach(@fos){
				if($mods_new[$_]){
					my $t=$mods_new[$_];
					$t=~s/[\(,\)]//g;
					$mods_new[$_] = '('.$t.',21)';	
				}
				else{
					$mods_new[$_] = '(21)';
				}
			}
			#Genero el peptido
			my $peptide_string = '';
			for(my $i=0; $i<$#aas+1; $i++){
				if($mods_new[$i]){
					$peptide_string = $peptide_string.$aas[$i].$mods_new[$i];
				}
				else{
					$peptide_string = $peptide_string.$aas[$i];
				}
			}
			push @combinaciones_peptidos,$peptide_string;
		}
	}
	elsif(@combinaciones_DH){
		#print "solo DH\t";
		foreach my $dhs(@combinaciones_DH){
			my @dhs=@$dhs;
			#Genero nuevo array de modificaciones para construir el peptido modificado
			my @mods_new;
			#Recupero las modificaciones no phospho o DH (chemtag,MetOx,...)
			for(my $i=0; $i<$#mods+1; $i++){
				if($mods[$i]){
					unless(($mods[$i]==21) or ($mods[$i]==23)){
						$mods_new[$i] = '('.$mods[$i].')';
					}
				}
			}
			foreach(@dhs){
				if($mods_new[$_]){
					my $t=$mods_new[$_];
					$t=~s/[\(,\)]//g;
					$mods_new[$_]='('.$t.',23)';
				}
				else{
					$mods_new[$_]='(23)';
				}
			}
			#Genero el peptido
			my $peptide_string='';
			for(my $i=0;$i<$#aas+1;$i++){
				if($mods_new[$i]){
					$peptide_string = $peptide_string.$aas[$i].$mods_new[$i];
				}
				else{
					$peptide_string = $peptide_string.$aas[$i];
				}
			}
			push @combinaciones_peptidos, $peptide_string;
		}
	}
	else{
		push @combinaciones_peptidos,$_[0];
	}
	return @combinaciones_peptidos;	
}

sub prever_fragmentos{
	# Module to predict a spectra using a sequence. It uses the 20 aa and allows phosphorilations on STY (+80:s,t,y),
	#water loss on ST (-18:u,v), and methionine oxidation (+18:m). It produces a hash containing the {label}->mass.
	#It calculates the y and b series for the charge state, 
	#ARGUMENTS: there are 3 args:
	#       - $_[0]-> sequence
	#       - $_[1]-> charge state (1,2 or 3)
	#       - $_[2]-> mz upper limit for fragment generation
	#RETURN: a hash with {label}->mass.
	my $peptide = $_[0];
	my $z = $_[1];
	my $limit_mz_sup = $_[2];
	my $limit_mz_inf = $_[3];
	my %masas_aa = (
	A=>71.03711,  V=>99.06841,  L=>113.08406, I=>113.08406, P=>97.05276,
	M=>131.0405,  F=>147.0684,  W=>186.07932, G=>57.02146,  S=>87.03203,
	T=>101.04768, C=>160.03069, Y=>163.06332, N=>114.04293, Q=>128.05858,
	D=>115.02694, E=>129.04259, H=>137.05891, K=>128.09496, R=>156.1011);
	#la C esta CAM!, la considero PTM fija
	
	my @aas;
	my @mods;
	my @masas_aas;

	while($peptide=~/\w/){
		my $aa;
		my $mod;
		if($peptide=~s/^(\w)\((\w+\,\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
		}
		elsif($peptide=~s/^(\w)\((\w+)\)//){
			$aa = $1;
			$mod = $2;
			push @aas, $aa;
			push @mods, $mod;
		}
		else{
			$peptide=~s/^(\w)//;
			$aa = $1;
			$mod = '';
			push @aas,$aa;
			push @mods,$mod;
		}
	}
	my %masas_mod = (21=>79.9799,         #fosfato
                     23=>-18.010565,      #dh
                     214=>144.102063,     #itraq
                     737=>229.162932,     #tmt
                     35=>15.994915);      #oxidacion

        for(my $i=0;$i<$#aas+1;$i++){
		if($mods[$i]){
			if($mods[$i]=~/\,/){
				my($m1,$m2)=split /\,/,$mods[$i];
				$masas_aas[$i] = $masas_mod{$m1} + $masas_mod{$m2} + $masas_aa{$aas[$i]};
			}
			else{
				$masas_aas[$i] = $masas_mod{$mods[$i]} + $masas_aa{$aas[$i]};
				#print $mods[$i],"\n";
			}
		}
		else{
			$masas_aas[$i] = $masas_aa{$aas[$i]};
		}
	}
	
	#La carga maxima permitida es 1 para z=1 o 2, 
	my $fragment_max_charge;
	if($z<3){
		$fragment_max_charge = 1;
	}
	else{
		$fragment_max_charge = $z-1;
	}

	my $acumulado = 0;
	my (@serie_y_1, @serie_b_1, @serie_y_2, @serie_b_2, @serie_y_3, @serie_b_3);

	for (my $i=0; $i<$#aas+1; $i++){
		$serie_b_1[$i] = $masas_aas[$i] + $acumulado+1;
		$acumulado += $masas_aas[$i];
		}
	my $M = $acumulado+19;								####ojo este 19!!!!!!
	for(my $i=0;$i<$#aas;$i++){
		$serie_y_1[$i] = $M-$serie_b_1[($#aas-$i)]+1;
		}
	$serie_y_1[$#aas] = $M-$serie_b_1[0]+1;
	if($_[1]>1){
		for (my $i=0; $i<$#aas+1; $i++){
			$serie_y_2[$i] = ($serie_y_1[$i]+1)/2;
			if (defined $serie_b_1[$i]){
				$serie_b_2[$i] = ($serie_b_1[$i]+1)/2;
				}
			}
		if($_[1]>2){
			for (my $i=0; $i<$#aas+1; $i++){
				$serie_y_3[$i] = ($serie_y_1[$i]+1)/3;
				$serie_b_3[$i] = ($serie_b_1[$i]+1)/3;
			}		
		}	    	
}

	my %fragmentos;
	for (my $i=0;$i<$#aas+1;$i++){
		my $j=$i+1;
		if ($i<($#aas)){															# El b(n) bo existe
			my $num = $j;
			if($num<10){$num = '0'.$num};
			$fragmentos{"b$num+1"} = sprintf("%.3f", $serie_b_1[$i]);
			if(($_[1]>1)and($fragment_max_charge>1)){					# Si carga>1
				$fragmentos{"b$num+2"} = sprintf("%.3f", $serie_b_2[$i]);
				if(($_[1]>2) and ($fragment_max_charge>2)){			# Si carga>2
					$fragmentos{"b$num+3"} = sprintf("%.3f", $serie_b_3[$i]);
				}
			}	
		}
		if(($i>0)){	
			my $num = $i;
			if($num<10){$num='0'.$num};												# El y(0) bo existe
			$fragmentos{"y$num+1"} = sprintf("%.3f", $serie_y_1[$i]);
			 if(($_[1]>1)and($fragment_max_charge>1)){				# Si carga>1
				$fragmentos{"y$num+2"} = sprintf("%.3f", $serie_y_2[$i]);
				if(($_[1]>2)and($fragment_max_charge>2)){				# Si carga>2
					$fragmentos{"y$num+3"} = sprintf("%.3f", $serie_y_3[$i]);	
				}	
			 }
		}	
	}
	#Solo se consideran los fragmentos entre limit_mz_sup e _inf
	foreach my $frag(keys %fragmentos){
		if(($fragmentos{$frag}>$limit_mz_sup)or($fragmentos{$frag}<$limit_mz_inf)){
			delete $fragmentos{$frag};
		}
	}
	return \%fragmentos;
}


sub matches_10niveles_series_yb_completas{
	my %fragmentos = %{$_[0]};
	my @masas = @{$_[1]};
	my @intensidades = @{$_[2]};
	my $fragment_error = $_[3];
	my $peptido = $_[4];
	
	my $num_fragmentos = 0;
	foreach(keys %fragmentos){
		$num_fragmentos++;
		#print $_,"\t",$fragmentos{$_},"\n";
	}
	
	my $num_picos_cada_100;
	my %datos_num_picos;
	my @rangos_masas;
	my ($rango_masas_inf,$rango_masas_sup) = (200, 2000);
	for(my $i=$rango_masas_inf; $i<$rango_masas_sup; $i=$i+100){
		push @rangos_masas, $i;
	}
	for($num_picos_cada_100=1;$num_picos_cada_100<11;$num_picos_cada_100++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias = 0;
		my $num_picos_dentro_intervalo = 0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior=$intervalo-100;
			for(my $i=0;$i<$#masas+1;$i++){
				if(($masas[$i]>$intervalo_inferior) and ($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]} = $intensidades[$i];
					$num_picos_dentro_intervalo++;
				}
			}
			my $num_tests_dentro_intervalo=0;
			my @masas_ordenadas_por_intensidad;
			
			if($num_picos_dentro_intervalo==0){#ningun pico experimental en el intervalo
										       #No hay picos en el intervalo\n";
				next;
			}
			elsif($num_picos_dentro_intervalo==1){ #un solo pico experimental en el intervalo
				foreach my $key(keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$key;
				}
			}
			else{#mas de 1 pico experimental en el intervalo
				foreach my $k (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$k;
				}
			}
			
			#	Aqui se realiza el matching entre masas teoricas y experimentales
			my $masas_testadas = 0;
			foreach my $masas_exp(@masas_ordenadas_por_intensidad){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp)<=$fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo,$masas_exp;
						push @masas_teor_match_intervalo,$fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas==$num_picos_cada_100){
					last;
				}
			}
		}

		my $probablity = ($num_picos_cada_100/100);
		my $binomial_acum_rango = binomial($num_fragmentos, $correspondencias, $probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score = -10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score = 0;
		}
		
		$datos_num_picos{$num_picos_cada_100}{masas_exp_match_intervalo} = \@masas_exp_match_intervalo;
		$datos_num_picos{$num_picos_cada_100}{masas_teor_match_intervalo} = \@masas_teor_match_intervalo;
		$datos_num_picos{$num_picos_cada_100}{tags_match_intervalo} = \@tags_match_intervalo;
		$datos_num_picos{$num_picos_cada_100}{score} = $score;
	}
	
	return \%datos_num_picos;
}

sub binomial{
	#Binomial distribution
	#($n,$k,$p)->$k successes in $n tries, given a probability of $p	
	my $n = $_[0];
	my $k = $_[1];
	my $p = $_[2];
	my $binomial_acumulada=0;
	for(my $i=$k; $i<$n+1; $i++){
		my $binomial = Math::NumberCruncher::Binomial($n,$i,$p);
		$binomial_acumulada = $binomial_acumulada+$binomial;
	}
	return $binomial_acumulada;
}

sub conversion_styuv{
	my $pep=$_[0];
	#cambio fosfo o DH simples
	$pep=~s/S\(21\)/s/g;
	$pep=~s/T\(21\)/t/g;
	$pep=~s/Y\(21\)/y/g;
	$pep=~s/S\(23\)/u/g;
	$pep=~s/T\(23\)/v/g;
	#cambio fosfo o DH combinados con chemtag en Nt
	$pep=~s/S\((214|737)\,21\)/s/;
	$pep=~s/T\((214|737)\,21\)/t/;
	$pep=~s/Y\((214|737)\,21\)/y/;
	$pep=~s/S\(21\,(214|737)\)/s/;
	$pep=~s/T\(21\,(214|737)\)/t/;
	$pep=~s/Y\(21\,(214|737)\)/y/;
	$pep=~s/S\((214|737)\,23\)/u/;
	$pep=~s/T\((214|737)\,23\)/v/;
	$pep=~s/S\(23\,(214|737)\)/u/;
	$pep=~s/T\(23\,(214|737)\)/v/;
	#Elimino todas las demas modificaciones
	$pep=~s/\(\d+\)//g;
	return $pep;
}

sub tratamiento_fragmentos{
	#argumentos-> $_[0]:%fragmentos,$_[1]:$peptido,$_[2]:primer aa intervalo, $_[3]: ultimo aa intervalo;
	########		Tratamiento fragmentos		######
	#Genero un nuevo hash para poner los fragmentos que estan entre el primero y ultimo aa que se pueda fosforilar
		
	my $peptido = $_[0];
	my %fragmentos = %{$_[1]};
	#print Dumper %fragmentos;
	#exit;
	my $primer_aa_fosforilable = $_[2];
	my $ultimo_aa_fosforilable = $_[3];
	#print "start-end, $_[2], $_[3]\n";
	my %fragmentos_filtrados;

	my $largo_pep = length $peptido;
	
	foreach my $frag (keys %fragmentos){
		my $frag2=(split /\+/,$frag)[0];
		if($frag2=~m/y/){
			$frag2=~s/y//;		#para las y, valen los fragmentos entre [largo-ultimo+1 , largo-primero]
			if(($frag2 >= ($largo_pep-$ultimo_aa_fosforilable+1))and
			   ($frag2 <= ($largo_pep-$primer_aa_fosforilable))){
				$fragmentos_filtrados{$frag} = $fragmentos{$frag};
			}
		}
		elsif($frag2=~m/b/){	#para las b, valen los fragmentos entre [primero , ultimo-1]
			$frag2=~s/b//;
			if(($frag2>=($primer_aa_fosforilable))and($frag2<=($ultimo_aa_fosforilable-1))){
				$fragmentos_filtrados{$frag}=$fragmentos{$frag};
			}
		}
	}
	%fragmentos = ();
	%fragmentos = %fragmentos_filtrados;	#Uso solo fragmentos en zona de interes		
	#######		fin tratamiento fragmentos
	return \%fragmentos;
}


sub matches_ascore{
	my %fragmentos = %{$_[0]};
	my @masas= @{$_[1]};
	my @intensidades = @{$_[2]};
	my $fragment_error = $_[3];
	my $peptido = $_[4];
	
	#print "ascore_frag_error, $_[3]\n";
	#print "peptido, $peptido\n";
	
	my $num_fragmentos = 0;
	foreach(keys %fragmentos){
		$num_fragmentos++;
	}
	
	my $num_picos_cada_100;
	my %datos_num_picos;
	my @rangos_masas;
	my ($rango_masas_inf, $rango_masas_sup) = (200, 2000);
	for(my $i=$rango_masas_inf; $i<$rango_masas_sup; $i=$i+100){
		push @rangos_masas, $i;
	}
	for($num_picos_cada_100=1; $num_picos_cada_100<11; $num_picos_cada_100++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;
		my $correspondencias = 0;
		my $num_picos_dentro_intervalo = 0;
		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior = $intervalo-100;
			for(my $i=0; $i<$#masas+1; $i++){
				if(($masas[$i]>$intervalo_inferior) and ($masas[$i]<=$intervalo)){
					$picos_intervalo{$masas[$i]} = $intensidades[$i];
					$num_picos_dentro_intervalo++;
				}
			}
			my $num_tests_dentro_intervalo=0;
			my @masas_ordenadas_por_intensidad;
			
			#ningun pico experimental en el intervalo
			if($num_picos_dentro_intervalo==0){
				next;
			}
			#un solo pico experimental en el intervalo
			elsif($num_picos_dentro_intervalo==1){ 
				foreach my $key(keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$key;
				}
			}
			#mas de 1 pico experimental en el intervalo
			else{	
				foreach my $k (sort { $picos_intervalo{$b} <=> $picos_intervalo {$a}} keys %picos_intervalo ){
					push @masas_ordenadas_por_intensidad,$k;
				}
			}
			
			#Matching entre masas teoricas y experimentales
			my $masas_testadas=0;
			foreach my $masas_exp(@masas_ordenadas_por_intensidad){
				$masas_testadas++;
				foreach(keys %fragmentos){
					if(abs($fragmentos{$_}-$masas_exp)<=$fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo,$masas_exp;
						push @masas_teor_match_intervalo,$fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas==$num_picos_cada_100){
					last;
				}
			}
		}
		#$num_fragmentos=> trials
		#$correspondencias=>  successes
		#$prob=> numero picos cada 100 / 100

		my $probablity = ($num_picos_cada_100/100);
		my $binomial_acum_rango = binomial($num_fragmentos, $correspondencias, $probablity);
		my $score;
		if($binomial_acum_rango>0){
			$score = -10*(log($binomial_acum_rango)/log(10));
		}
		else{
			$score = 0;
		}
		#print "score, $score\n";
		
		$datos_num_picos{$num_picos_cada_100}{masas_exp_match_intervalo} = \@masas_exp_match_intervalo;
		$datos_num_picos{$num_picos_cada_100}{masas_teor_match_intervalo} = \@masas_teor_match_intervalo;
		$datos_num_picos{$num_picos_cada_100}{tags_match_intervalo} = \@tags_match_intervalo;
		$datos_num_picos{$num_picos_cada_100}{score} = $score;
	}
	
	return \%datos_num_picos;

}


return 1;

__END__

=head1 NAME

=head1 SYNOPSIS

=head1 DESCRIPTION

=head1 LIMITATIONS

=head1 SEE ALSO

=cut



